/*		
	
	This header-file is filled with functions from solutions to previous assignments:
	-filter.h & filter.cpp
	-graymorphology.h
	-MoreFiltering.h
	This source was created by supervisors of the image processing courses.
	As these are exemplarely solutions, only changed lines of code are commented.
	Some unneeded functionalities were removed without comment.

*/

#ifndef __IP_H__
#define __IP_H__

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
#include <cassert>
#include <list>
#include <cmath>

#include "Image.H"
#include "IO.H"
#include "Util.H"

//
//   code of graymorphology.h:
//

void discStructureElement(pip::Image<unsigned char> &imOut, int rad){

     pip::Dimension dim(2*rad+1,2*rad+1); // dimensions of SE --> +1 for center
     imOut.resize(dim);
     int rad2 = rad*rad;
     int val;
     for (int y=0; y<=rad; ++y){
        for (int x=0; x<=rad; ++x){
          val = 0;

          if( (x*x+y*y) < rad2+0.5)
              val = 1;

          imOut(rad+x,rad+y) = val;
          imOut(rad-x,rad+y) = val;
          imOut(rad-x,rad-y) = val;
          imOut(rad+x,rad-y) = val;
      }
	}
}

template <class Tin>
void grayDilate(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut, const pip::Image<unsigned char> &se){
        // make sure we have the correct output size
        if (imOut.dimension()!= imIn.dimension()) {
                imOut.resize(imIn.dimension());
        }

        int iSEWidth  = se.dim(0);
        int iSECenter = iSEWidth/2;
        int iWidth    = imIn.dim(0);
        int iHeight   = imIn.dim(1);

        int iX, iY;
        // over all image pixels
        for (int y=0; y<iHeight; ++y)
        for (int x=0; x<iWidth; ++x) {
                imOut(x,y) = imIn(x,y);
                // over all se pixels
                for (int m=0; m<iSEWidth; ++m)
                for (int n=0; n<iSEWidth; ++n) {
                        if (se(n,m) > 0) {
                                iX = x+(m-iSECenter);
                                iY = y+(n-iSECenter);
                                if ( iX >= 0 && iX < iWidth && iY >= 0 && iY < iHeight)
                                imOut(x, y) = max(imOut(x,y),imIn(iX,iY));
                        }
                }
        }
}

template <class Tin>
void grayErode(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut, const pip::Image<unsigned char> &se){
        // make sure we have the correct output size
        if (imOut.dimension()!= imIn.dimension()) {
                imOut.resize(imIn.dimension());
        }

        int iSEWidth  = se.dim(0);
        int iSECenter = iSEWidth/2;
        int iWidth    = imIn.dim(0);
        int iHeight   = imIn.dim(1);

        int iX, iY;
        // over all image pixels
        for (int y=0; y<iHeight; ++y)
        for (int x=0; x<iWidth; ++x) {
                imOut(x,y) = imIn(x,y);
                // over all se pixels
                for (int m=0; m<iSEWidth; ++m)
                for (int n=0; n<iSEWidth; ++n) {
                        if (se(n,m) > 0) {
                                iX = x+(m-iSECenter);
                                iY = y+(n-iSECenter);
                                if ( iX >= 0 && iX < iWidth && iY >= 0 && iY < iHeight)
                                imOut(x, y) = min(imOut(x,y),imIn(iX,iY));
                        }
                }
        }
}

template <class Tin>
void grayOpening(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut, const pip::Image<unsigned char> &se){
        grayErode(imIn, imOut,se);
        pip::Image<Tin> tmpCopy(imOut);
        grayDilate(tmpCopy, imOut,se);
}

template <class Tin>
void grayClosing(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut, const pip::Image<unsigned char> &se){
        grayDilate(imIn, imOut,se);
        pip::Image<Tin> tmpCopy(imOut);
        grayErode(tmpCopy, imOut,se);
}


//
//   end of code of graymorphology.h
//


//
//   code of filtering.h:
//

void convertFloatToUnsignedChar(const pip::Image<float> &input, pip::Image<unsigned char> &output);
void convertUnsignedCharToFloat(const pip::Image<unsigned char> &input, pip::Image<float> &output);

float linearInterpolation(const pip::Image<float> &image, const float &x, const float &y);

void filter(const pip::Image<float> &input, pip::Image<float> &output, const pip::Image<float> &mask, const pip::BorderInfo &bi = pip::border_mirror);
void filterSeparable(const pip::Image<float> &in, pip::Image<float> &out, const std::vector<std::vector<float> > &kernels, const pip::BorderInfo &bi = pip::border_mirror);

float g0(const int &i, const float &sigma);
float g1(const int &i, const float &sigma);
float g2(const int &i, const float &sigma);
float g3(const int &i, const float &sigma);

void normalize(std::vector<float> &kernel);

void gauss(const pip::Image<float> &in, pip::Image<float> &out, const float &sigma, const int &dx = 0, const int &dy = 0, const float &cutoff = 5.0, const pip::BorderInfo &bi = pip::border_mirror);

// Return minimum and maximum of the input image
template <class T>
void getMinMaxImage(const pip::Image<T> &in, T &min, T &max)
{
  min = numeric_limits<T>::max();
  max = numeric_limits<T>::min();

  for (int z = 0; z < in.dim(2); ++z) {
    for (int y = 0; y < in.dim(1); ++y) {
      for (int x = 0; x < in.dim(0); ++x) {
        T value(in(x, y, z));

        if (value > max) {
          max = value;
        }

        if (value < min) {
          min = value;
        }
      }
    }
  }
}

//
//   end of code of filter.h
//


//
//   code of MoreFiltering.h:
//

/// MEDIAN FILTERING \\\

/**
 * Performs median filtering on an image, using mirrored borders.
 *
 * @param image Image to filter
 * @param n Size of the mask (kernel)
 * @param squareMask Whether the mask is a square shape or a plus shape
 */
void medianFilter(Image<unsigned char> &image, const int &n, const bool &squareMask = true)
{
  const int kernelCentre(n / 2);
	Image<unsigned char> tempImage(image.bordered(kernelCentre, border_mirror));

  // Vector used to collect pixels over which median will be computed
  vector<unsigned char> pixels;

  for (int y(0); y < image.dim(1); ++y) {
    for (int x(0); x < image.dim(0); ++x) {
      for (int dy(-kernelCentre); dy <= kernelCentre; ++dy) {
        for (int dx(-kernelCentre); dx <= kernelCentre; ++dx) {
          // If we use a square mask, all pixels need to be processed.
          // If we use a plus mask, only pixels in the centre need to be processed.
          if (squareMask || dy == 0 || dx == 0) {
            pixels.push_back(tempImage(x + dx, y + dy));
          }
        }
      }

      // Iterator pointing to the median element of the vector.
	    vector<unsigned char>::iterator medianIterator(pixels.begin() + pixels.size() / 2);

      nth_element(pixels.begin(), medianIterator, pixels.end());
	    image(x, y) = *medianIterator;

      // Clear the vector, but keep its capacity
      pixels.erase(pixels.begin(), pixels.begin() + pixels.size());
    }
  }
}

/**
 * Performs fast median filtering on an image, using mirrored borders. The implementation is based
 * on the algorithm in Sonka, page 129-130. Adapted by Hugo Kuijf for extra speed.
 * It works for square kernels only, for plus-shaped kernels the above medianFilter method is called.
 * The same principle can be applied to plus-kernels, although the gain would not be that big as a
 * larger part of the kernel changes during the shift.
 *
 * Timings:
 * medianFilter as above        : 0.506 seconds
 * medianFilterFast as in Sonka : 0.032 seconds
 * medianFilterFast as below    : 0.027 seconds
 *
 * @param image Image to filter
 * @param n Size of the mask (kernel)
 * @param squareMask Whether the mask is a square shape or a plus shape
 */
void medianFilterFast(Image <unsigned char> &image, const int &n, const bool &squareMask = true)
{
  if (!squareMask) {
    return medianFilter(image, n, squareMask);
  }

  // Initialize variables
  int dx, dy, x, y, kernelCentre(n / 2), nm, initialNm, t(n * n / 2);
  Image<unsigned char> tempImage(image.bordered(kernelCentre, border_mirror));
  unsigned char m, initialM, pg;

  // Create initial histogram at (0,0) and compute the median
  int histogram[256], initialHistogram[256] = {0};

  //   Fill histogram
  for (dy = -kernelCentre; dy <= kernelCentre; ++dy) {
    for (dx = -kernelCentre; dx <= kernelCentre; ++dx) {
      initialHistogram[tempImage(dx, dy)]++;
    }
  }

  //   Compute initial mean (m) and number of pixels before the mean (nm)
  for (initialM = 0, initialNm = initialHistogram[initialM]; initialNm <= t; initialNm += initialHistogram[++initialM]);

  // Process whole image
  for (y = 0; y < image.dim(1); ++y) {
    // Copy the initial histogram, mean and nm for this row to the current histogram, mean and nm
    memcpy(histogram, initialHistogram, 256 * sizeof(int));
    m = initialM; nm = initialNm;

    // Shift histogram down. Only not at the first row, then we already have the correct histogram.
    if (y != 0) {
      // Create new initial histogram for this row, based on the initial histogram of the previous row
      for (dx = -kernelCentre; dx <= kernelCentre; ++dx) {
        // Clear top column
        pg = tempImage(dx, (y - 1) - kernelCentre);
        histogram[pg]--;
        if (pg <= m) --nm;

        // Fill bottom column
        pg = tempImage(dx, y + kernelCentre);
        histogram[pg]++;
        if (pg <= m) ++nm;
      }

      // Find the median. The median is at the location where the cumulative histogram (stored in nm) is equal to or larger than t.
      for (; nm <= t; ++m, nm += histogram[m]);
      for (; nm - histogram[m] > t; nm -= histogram[m], --m);

      // Copy the current histogram for this row, so the next row can use it as the initial histogram
      memcpy(initialHistogram, histogram, 256 * sizeof(int));
      initialM = m; initialNm = nm;
    }

    // Set the median
    image(0, y) = m;

    for (x = 1; x < image.dim(0); ++x) {
      for (dy = -kernelCentre; dy <= kernelCentre; ++dy) {
        // Clear left column
        pg = tempImage((x - 1) - kernelCentre, y + dy);
        histogram[pg]--;
        if (pg <= m) --nm;

        // Fill right column
        pg = tempImage(x + kernelCentre, y + dy);
        histogram[pg]++;
        if (pg <= m) ++nm;
      }

      // Find the median. The median is at the location where the cumulative histogram (stored in nm) is equal to t.
      for (; nm <= t; ++m, nm += histogram[m]);
      for (; nm - histogram[m] > t; nm -= histogram[m], --m);

      // Set the median
      image(x, y) = m;
    }
  }
}

/// SHARPENING \\\

/**
 * Initializes the Laplacian kernel
 *
 * @param kernel The Laplacian kernel
 * @param scale
 */
void initLaplacianKernel(Image<float> &kernel, const float &scale)
{
  kernel.resize(Dimension(3, 3));

  kernel(0,0) = -scale;
  kernel(1,0) = -scale;
  kernel(2,0) = -scale;

  kernel(0,1) = -scale;
  kernel(1,1) = 1 + 8 * scale;
  kernel(2,1) = -scale;

  kernel(0,2) = -scale;
  kernel(1,2) = -scale;
  kernel(2,2) = -scale;
}

/**
 * Sharpens the image by subtracting a Laplacian with factor C.
 * See also 5.35 in Sonka
 *
 * @param image Image to sharpen
 * @param C Factor
 */
void sharpen(Image<float> &image, const float &C)
{
  // Initialize the filter
  Image<float> kernel;
  initLaplacianKernel(kernel, C);

  Image<float> out;
  filter(image, out, kernel);

  image = out;
}

/// UNSHARP MASKING \\\

/**
 * Performs unsharp masking on an image.
 * 
 * @param image Image to filter
 * @param sigma
 * @param a
 */
void unsharpMasking(Image<float> &image, const float &sigma, const float &a)
{
  Image<float> gaussImage;
  gauss(image, gaussImage, sigma, 0, 0, 5.0);

  for (int y = 0; y < image.dim(1); ++y) {
    for (int x = 0; x < image.dim(0); ++x) {
      image(x, y) = (1.0f - a) * gaussImage(x, y) + a * image(x, y);
    }
  }
}

/// DENOISE \\\

/**
 * Generates a random number in the range [min, max]
 *
 * @param min Lower range value
 * @param max Upper range value
 * @return Random number
 */
float randomNumber(const float &min, const float &max)
{
  return min + ((max - min) * rand() / RAND_MAX);
}

/**
 * Generate a random number from a Gaussian distribution N(mean, sigma)
 *
 * @param mean
 * @param sigma
 * @return Random number
 */
float randomGaussian(const float &mean, const float &sigma)
{
  float U1, U2, V1, V2, S;

  do {
    U1 = randomNumber(0, 1);
    U2 = randomNumber(0, 1);
    V1 = 2 * U1 - 1;
    V2 = 2 * U2 - 1;
    S = V1 * V1 + V2 * V2;
  } while (S >= 1);

  float X = sqrt(-2 * log(S) / S) * V1;

  return (mean + sigma * X);
}

/**
 * Degrade an image with additive Gaussian noise.
 *
 * @param image Image to degrade
 * @param sigma
 */
void degradeGauss(Image<unsigned char> &image, const float &sigma)
{
  for (int y = 0; y < image.dim(1); ++y) {
    for (int x = 0; x < image.dim(0); ++x) {
      float randGauss = randomGaussian(0, sigma);

      int sign = randGauss < 0 ? -1 : 1,
        greyVal = (int)floor(abs(randGauss) + 0.5),
        noisyValue = image(x, y) + sign * greyVal;

      if (noisyValue < 0) {
        image(x, y) = 0;
      } else if (noisyValue > 255) {
        image(x, y) = 255;
      } else {
        image(x, y) = (unsigned char) noisyValue;
      }
    }
  }
}

/**
 * Degrade an image with impulse noise. Fraction denotes the
 * fraction of pixels affected by corruption.
 *
 * @param image Image to degrade
 * @param fraction
 */
void degradeImpulse(Image<unsigned char> &image, const float &fraction)
{
  for (int y = 0; y < image.dim(1); ++y) {
    for (int x = 0; x < image.dim(0); ++x) {
      float randUniform = randomNumber(-1, 1);

      if (abs(randUniform) < fraction) {
        image(x, y) = randUniform < 0 ? 0 : 255;
      }
    }
  }
}

/**
 * Computes the average squared difference between two input images.
 *
 * @param input First input image
 * @param secondInput Second input image
 * @return Average squared difference
 */
float averageSquaredDifference(const Image<unsigned char> &input, const Image<unsigned char> &secondInput)
{
  float diff(0);
  for (int y = 0; y < input.dim(1); ++y) {
    for (int x = 0; x < input.dim(0); ++x) {
      diff += (input(x,y) - secondInput(x,y)) * (input(x,y) - secondInput(x,y));
    }
  }

  return (diff / (input.dim(0) * input.dim(1)));
}

/**
 * Degrades an image with noise (additive Gaussian or Impulse noise),
 * and restores it again (by Gaussian blurring or median filtering).
 *
 * @param image Input image
 * @param degradedImage Output degraded image
 * @param degradedFilteredImage The 'restored' degraded output image
 * @param s1 Parameter for degrading (sigma for Gaussian noise or fraction for impulse noise)
 * @param s2 Parameter for restoring (sigma for blurring or n for median filtering)
 * @param filterMethod True for Gaussian blurring, false for median filtering
 */
void denoise(const Image<unsigned char> &image,
             Image<unsigned char> &degradedImage, 
             Image<unsigned char> &degradedFilteredImage, 
             const float &s1, 
             const float &s2, 
             const bool &filterMethod,
             void (*additiveNoise)(Image<unsigned char>&, const float &))
{
  // Make sure to check if filterMethod has an allowed value
  // Make sure that s1 is in a valid range, depending on the noise method used
  if(additiveNoise == degradeGauss && s1 <= 0) {
    cerr << "Incorrect value for s1, sigma must be > 0" << endl;
    return;
  } else if (additiveNoise == degradeImpulse && (s1 < 0 || s1 > 1)) {
    cerr << "Incorrect value for s1, sigma must be > 0 and < 1" << endl;
    return;
  }

  // Make sure s2 has valid value, it should be an odd number for median filtering
  if(!filterMethod && (int)s2 % 2 == 0) { // median filtering with even kernel size
    cerr << "Median filtering only allowed with odd kernel size" << endl;
    return;
  }

  additiveNoise(degradedImage, s1);

  if(filterMethod) { // Gaussian averaging
    Image<float> gaussImage, tempImage;

    convertUnsignedCharToFloat(degradedImage, tempImage);
    gauss(tempImage, gaussImage, s2);
    convertFloatToUnsignedChar(gaussImage, degradedFilteredImage);
  } else { // Median filtering
    degradedFilteredImage = degradedImage;
    medianFilter(degradedFilteredImage, (int)s2);
  }
}

/// CANNY EDGE DETECTION \\\

/**
 * Gradient magnitude ...
 *
 * @param in Input image
 * @param dx Output of X-derivative
 * @param dy Output of Y-derivative
 * @param gm Output of Gradient Magnitude
 * @param sigma
 */
void gradientMagnitude(const Image<unsigned char> &in, Image<float> &dx, Image<float> &dy, Image<float> &gm, const float &sigma)
{
  Image<float> image;
  gm.resize(in.dimension());

  convertUnsignedCharToFloat(in, image);

  gauss(image, dx, sigma, 1);
  gauss(image, dy, sigma, 0, 1);

  for(int y = 0; y < gm.dim(1); ++y) {
    for(int x = 0; x < gm.dim(0); ++x) {
      gm(x, y) = sqrt(dx(x, y) * dx(x, y) + dy(x, y) * dy(x, y));
    }
  }
}

/**
 * Canny edge detector ...
 *
 * @param in Input image
 * @param out Output of gradient magnitude for pixels that survive non-maximal suppression
 * @param sigma Scale for the gradient magnitude
 */
void cannyEdgeDetector(const Image<unsigned char> &in, Image<float> &out, const float &sigma)
{
  out.resize(in.dimension());
  out.fill(0);

  Image<float> dxImage, dyImage;
  Image<float> gmImage(in.dimension());

  // Compute gradient magnitude and edge direction
  gradientMagnitude(in, dxImage, dyImage, gmImage, sigma);

  // Apply non-maximum suppression in direction of the edge
  for(int y = 0; y < in.dim(1); ++y) {
    for(int x = 0; x < in.dim(0); ++x) {
      float gMag(gmImage(x, y));

      //if(gMag > 1.0E-10) { //if smaller, set to 0
        // calculate normal to edge according to formula 4.63 in Sonka
        // normal = derivative/gradientmagnitude
        float dx(dxImage(x, y) / gMag),
          dy(dyImage(x, y) / gMag);

        // do non-maximal suppression --> get only the strongest edge
        // Get gradientmagnitude interpolation value in positive and negative direction normal to edge
        // to check if these edges are stronger!
        float gMag1(linearInterpolation(gmImage, x + dx, y + dy)),
          gMag2(linearInterpolation(gmImage, x - dx, y - dy));

        if(gMag1 < gMag && gMag2 < gMag) { // x, y is maximum
          out(x, y) = gMag;
       // }
      } else {
        out(x, y) = 0;
      }
    }
  }
}

/**
 * Number of edge voxels ...
 *
 * @param filteredImage ...
 * @return ...
 */
int nrOfEdgeVoxels(const Image<float> &filteredImage)
{
  int sum(0);

  for(int y = 0; y < filteredImage.dim(1); ++y) {
    for(int x = 0; x < filteredImage.dim(0); ++x) {
      if(filteredImage(x, y) != 0) {
        ++sum;
      }
    }
  }

  return sum;
}

//
//   end of code of MoreFiltering.h
//

/**
 * Convert a Float image to an Unsigned Char image and scales the values between 0 and 255.
 * A constant float image is filled with value 128.
 *
 * @param input Constant reference to the input Float image
 * @param output Reference to the output Unsigned Char image
 */
void convertFloatToUnsignedChar(const Image<float> &input, Image<unsigned char> &output)
{
  float min, max;
  getMinMaxImage(input, min, max);

  output.resize(input.dimension());

  if (min == max) {
    output.fill(128);
    return;
  }
  
  float scale(255.0f / (max - min));

  for (int z = 0; z < input.dim(2); ++z) {
    for (int y = 0; y < input.dim(1); ++y) {
      for (int x = 0; x < input.dim(0); ++x) {
        output(x, y, z) = (unsigned char)(scale * (input(x, y, z) - min));
      }
    }
  }
}

/**
 * Converts an Unsigned Char image to a Float image.
 *
 * @param input Constant reference to the input Unsigned Char image
 * @param output Reference to the output Float image
 */
void convertUnsignedCharToFloat(const Image<unsigned char> &input, Image<float> &output)
{
  output.resize(input.dimension());

  for (int z = 0; z < input.dim(2); ++z) {
    for (int y = 0; y < input.dim(1); ++y) {
      for (int x = 0; x < input.dim(0); ++x) {
        output(x, y, z) = (float)input(x, y, z);
      }
    }
  }
}

/**
 * Returns the value at (x, y) using linear interpolation. See also
 * Sonke, page 122 and 123, equation 5.22
 *
 * @param image Input image
 * @param x
 * @param y
 * @return Interpolated value at (x, y)
 */
float linearInterpolation(const Image<float> &image, const float &x, const float &y)
{
  int l((int)x), 
    k((int)y);

  // Check if coordinates are within the image
  // Note that the image.dim(*) - 1 is necessary to prevent the interpolation from going outside the image
  if (l < 0 || k < 0 || l >= image.dim(0) - 1 || k >= image.dim(1) - 1) {
    return 0.0;
  } else {
    float a(x - l),
      b(y - k);

    return (float)((1.0 - a) * (1.0 - b) * image(l, k)
      + a * (1.0 - b) * image(l + 1, k)
      + b * (1.0 - a) * image(l, k + 1)
      + a * b * image(l + 1, k + 1));
  }
}

/**
 * Convolves the input image with a kernel and stores the result in the output.
 *
 * @param input Constant reference to the input image
 * @param output Reference to the output image
 * @param mask Constant reference to the kernel
 */
void filter(const Image<float> &input, Image<float> &output, const Image<float> &mask, const BorderInfo &bi)
{
  // Only 2D n x n kernels, with n is odd.
  if (mask.dim(0) != mask.dim(1) || mask.dim(0) % 2 == 0 || mask.nrDims() != 2) {
    cerr << "Given kernel is incorrect. It is either not squared, not odd-sized, or not 2 dimensional." << endl;
    return;
  }

  int centre((mask.dim(0) - 1) / 2);
  Image<float> temp = input.bordered(centre, bi);

  output.resize(input.dimension());
  output.fill(0);

  // Loop over all pixels
  for (int y = 0; y < input.dim(1); ++y) {
    for (int x = 0; x < input.dim(0); ++x) {
      // Loop over all kernel values of mask
      for (int dy = -centre; dy <= centre; ++dy) {
        for (int dx = -centre; dx <= centre; ++dx) {
          output(x, y) += temp(x + dx, y + dy) * mask(centre - dx, centre - dy);
        }
      }
    }
  }
}

/**
 * Performs a separable convolution. The first dimension is convolved with the first
 * kernel in kernels, and so on.
 *
 * @param in Input image
 * @param out Output image
 * @param kernels List of kernels
 * @param bi Border information
 */
void filterSeparable(const Image<float> &in, Image<float> &out, const vector<vector<float> > &kernels, const BorderInfo &bi)
{
  if (in.nrDims() != kernels.size()) {
    cout << in.nrDims() <<"nrDim\n";
	cout << kernels.size() <<"kernelsize\n";
    cerr << "Incorrect dimensions." << endl;
    return;
  }

  // Find out what the largest kernel half size is,
  // we need a border of this size
  int n(0);
  for (int i = 0; i < (int)kernels.size(); ++i) {
    if ((int)kernels[i].size() / 2 > n) {
      n = kernels[i].size() / 2;
    }
  }

  Image<float> temp1(in.bordered(n, bi));
  Image<float> temp2(in.dimension(), n, bi);
  temp2.fill(0);
  out.resize(in.dimension());
  out.fill(0);

  // First along x-direction,
  // store result in temp2
  vector<float> kernel = kernels[0];
  n = kernel.size() / 2;
  for (int y = 0; y < in.dim(1); ++y) {
    for (int x = 0; x < in.dim(0); ++x) {
      for (int i = -n; i <= n; ++i) {
        temp2(x, y) += temp1(x - i, y) * kernel[i + n];
      }
    }
  }

  // Now along y-direction,
  // store result in out
  kernel = kernels[1];
  n = kernel.size() / 2;
  for (int y = 0; y < in.dim(1); ++y) {
    for (int x = 0; x < in.dim(0); ++x) {
      for (int i = -n; i <= n; ++i) {
        out(x, y) += temp2(x, y - i) * kernel[i + n];
      }
    }
  }
}

/**
 * Samples the unnormalized Gaussian for given sigma at integer values.
 *
 * @param i
 * @param sigma
 * @return Unnormalized Gaussian
 */
float g0(const int &i, const float &sigma)
{
  return (float)(exp(-(float)(i * i) / (2.0 * sigma * sigma)));
}

/**
 * Samples the unnormalized first derivative of the Gaussian for given sigma at integer values.
 *
 * @param i
 * @param sigma
 * @return Unnormalized first derivative of the Gaussian
 */
float g1(const int &i, const float &sigma)
{
  return (float)(-(float)i / (sigma * sigma) 
    * exp(-(float)(i * i) / (2.0 * sigma * sigma)));
}

/**
 * Samples the unnormalized second derivative of the Gaussian for given sigma at integer values.
 *
 * @param i
 * @param sigma
 * @return Unnormalized second derivative of the Gaussian
 */
float g2(const int &i, const float &sigma)
{
  return (float)(1.0 / (sigma * sigma) * ((float)(i * i) / (sigma * sigma) - 1)
    * exp(-(float)(i * i) / (2.0 * sigma * sigma)));
}

/**
 * Samples the unnormalized third derivative of the Gaussian for given sigma at integer values.
 *
 * @param i
 * @param sigma
 * @return Unnormalized third derivative of the Gaussian
 */
float g3(const int &i, const float &sigma)
{
  return (float)(1.0 / (sigma * sigma * sigma * sigma) * ((3.0 * (float)(i) - (float)(i * i * i) / (sigma * sigma)))
    * exp(-(float)(i * i) / (2.0 * sigma * sigma)));
}

/**
 * Normalizes the vector to unit absolute sum.
 *
 * @param kernel Kernel to normalize
 */
void normalize(vector<float> &kernel)
{
  float f(0.0);
  for (unsigned int i = 0; i < kernel.size(); ++i) {
    f += fabs(kernel[i]);
  }
  
  if (f > 0) {
    for (unsigned int i = 0; i < kernel.size(); ++i) {
      kernel[i] /= f;
    }
  }
}

/**
 * Computes a Gaussian derivative of scale sigma and order dx, dy of in, result stored in out.
 *
 * @param in Input image
 * @param out Output image
 * @param sigma
 * @param dx
 * @param dy
 * @param cutoff
 * @param bi Border information
 */
void gauss(const Image<float> &in, Image<float> &out, const float &sigma, const int &dx, const int &dy, const float &cutoff, const BorderInfo &bi)
{
  //determine length of kernels based on sigma and cutoff
  int n((int)(sigma * cutoff));

  //fill kernels. If dx or dy are not 0,1,2 or 3 give an error and return
  if (dx < 0 || dx > 3 || dy < 0 || dy > 3) {
    return;
  }

  vector<vector<float> > kernels(2);
  for (int i = -n; i <= n; ++i) {
    switch(dx) {
      case 0:
        kernels[0].push_back(g0(i, sigma));
        break;
      case 1:
        kernels[0].push_back(g1(i, sigma));
        break;
      case 2:
        kernels[0].push_back(g2(i, sigma));
        break;
      case 3:
        kernels[0].push_back(g3(i, sigma));
        break;
      default:
        break;
    }

    switch(dy) {
      case 0:
        kernels[1].push_back(g0(i, sigma));
        break;
      case 1:
        kernels[1].push_back(g1(i, sigma));
        break;
      case 2:
        kernels[1].push_back(g2(i, sigma));
        break;
      case 3:
        kernels[1].push_back(g3(i, sigma));
        break;
      default:
        break;
    }
  }

  //normalize the kernels
  normalize(kernels[0]);
  normalize(kernels[1]);

  //call filterseparable
  filterSeparable(in, out, kernels, bi);
}




#endif